// 处理边界cell，主要添加重心点进行细分。
//
#include<tool/libcell_tools_algorithm.h>


#include<Mesh/libcell_iterator.h>

#define quote lib_cell_quote

static template_v* libcell_tools_xifen_cell(
    Mesh* m, template_c* c )
{
    if(c==NULL||m==NULL){return NULL;}

    int dim=c->vertices[0]->point_size;

    double *p=(double*)malloc(sizeof(double)*dim);
    memset(p,0,sizeof(double)*dim);

    for(int i=0;i<c->vertices_size;i++)
    {
        for(int j=0;j<dim;j++)
        {
            p[j]+=c->vertices[i]->point[j];
        }
    }
    if(c->vertices_size>0)
    {
        for(int j=0;j<dim;j++)
        {
            p[j]/=(double)(c->vertices_size);
        }
    }

    Node* nodehfv=NULL;
    int vs_size=c->halffaces[0]->vertices_size;
    // template_hf**hfs=(template_hf**)malloc(sizeof(template_hf*)*c->halffaces_size);
    // memset(hfs,0,sizeof(template_hf*)*c->halffaces_size);
    for(int i=0;i<c->halffaces_size;i++)
    {
        template_hf* hf=c->halffaces[i];

        template_v** vs=(template_v**)malloc(sizeof(template_v*)*(hf->vertices_size+1) );
        memmove(vs,hf->vertices,sizeof(template_v*)*hf->vertices_size);
        nodehfv=node_overlying(nodehfv,vs);
    }

    m->delete_cell(m,*c,true);

    template_v* re=m->create_vertexv(m,p,dim);
    free(p);
    //int i=0;
    for(Node*nit=nodehfv;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v** vs=(template_v**)(nit->value);

        vs[vs_size]=re;

        m->create_cellv(m,vs,vs_size+1);
        free(vs);
        // template_f* f=m->create_facev(m,vs,vs_size);
        // template_hf* hf=m->create_halfface(m,f,vs,vs_size);

        // hfs[i]=hf;
        // i++;
    }

    free_node(nodehfv);


    // free(hfs);
    return re;
}

static int libell_tools_is_xifen_legal(Mesh* m,
    template_c*c)
{
    int num=0;
    for(int i=0;i<c->halffaces_size;i++)
    {
        template_hf* hf=c->halffaces[i];
        if(m->face_is_boundary(m,*(hf->face) ))
        {
            num++;
        }
        if(num>=2){return 1;}
    }
    // if(num>=c->halffaces_size-1)
    // {
    //     return 1;
    // }
    return 0;
}

void libcell_tools_xifen_mesh_boundary(Mesh* m)
{
    Node* nodec=NULL;
    for(auto cit= m->c_begin(m);cit!=m->c_end(m);cit++)
    {
        if(libell_tools_is_xifen_legal(m,quote(cit)))
        {
            nodec=node_overlying(nodec,quote(cit));
        }
    }

    for(Node* nit=nodec;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        libcell_tools_xifen_cell(m,c);
    }


    free_node(nodec);
}

static inline int is_edge_legacy_to_collapse(Mesh* m,template_f* f)
{
    // printf("begin id edge leagacy \n");
    LB_Vector_voidptr vec=Mesh_vv_beginn(m, *(f->vertices[0]));
    // LB_Vector_voidptr* vecp1=&vec;
    for(int i=1;i<f->vertices_size;i++)
    {
        if(vec.size==0){break;}
        // vec.intersection();
        LB_Vector_voidptr vec2=Mesh_vv_beginn(m, *(f->vertices[i]));
        LB_Vector_voidptr vec1= vec.intersection(&vec,&vec2);
        vec.clear(&vec);
        vec2.clear(&vec2);
        vec=vec1;

    }
    int size= vec.size;
    // free_node(n1);
    vec.clear(&vec);
    return size;
}
//
// 把v1 合并到v2
// 要求必须是三角形
//
static template_v* merging_two_vertices(Mesh* m,template_v* v1,template_v* v2)
{
    if(v1==v2){return NULL;}
    if(m==NULL){return v2;}
    template_v* vs[2]={v1,v2};
    template_f* f=m->get_facev(m,vs,2);
    // printf(" v1 v2:%d %d\n",v1->id,v2->id);
    if(f==NULL||is_edge_legacy_to_collapse(m,f)>2)
    {
        // printf("f is null\n");
        return NULL;
    }

    Int_RB_Tree* tree_c=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_c);
    // printf("begin vc\n");
    // for(Node* nit=v1->cells;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto vcit=m->vc_begin(m,*v1);vcit!=m->vc_end(m,*v1);vcit++)
    {

        template_c*c=quote(vcit);
        // template_c* c=(template_c*)(nit->value);
        tree_c->insert(tree_c,c->id,c);
    }
    // printf("end vc\n");

    Int_RB_Tree* tree_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_hf);

    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            template_hf* hf=quote(chfit);
            if(tree_hf->find(tree_hf,hf->face->id)==NULL )
            {tree_hf->insert(tree_hf,hf->face->id,hf);}
            else
            {tree_hf->erase(tree_hf,hf->face->id);}
        }
    }

    // for(Node* nit=v1->faces;nit!=NULL;nit=(Node*)(nit->Next))
    for(auto vfit=m->vf_begin(m,*v1);vfit!=m->vf_end(m,*v1);vfit++)
    {
        template_f*f=quote(vfit);
        // template_f* f=(template_f*)(nit->value);
        tree_hf->erase(tree_hf,f->id);
    }

    // printf("hrere \n");
    // int cols=f->vertices[0]->point_size;
    int rows=f->vertices_size+1;
    // template_v*re=m->create_vertexv(m,p,cols);
    Node* node_hf=NULL;
    for(auto it=tree_hf->begin(tree_hf);it.it!=NULL;it++)
    {
        template_hf* hf=(template_hf*)(it.second);
        template_v** vs_value=(template_v**)malloc(sizeof(template_v*)*(hf->vertices_size+1) );
        memmove(vs_value,hf->vertices,sizeof(template_v*)*hf->vertices_size);
        vs_value[hf->vertices_size]=v2;
        node_hf=node_overlying(node_hf,vs_value);
        //it.it->data=vs_value;
    }

    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c*c=(template_c*)(it.second);
        m->delete_cell(m,*c,true);
    }
    // printf("hrere1 \n");

    for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v**vs_value=(template_v**)(nit->value);
        if(vs_value[0]==v2||vs_value[1]==v2)
        {
            free(vs_value);
            continue;
        }


        m->create_cellv(m,vs_value,rows);
        free(vs_value);
    }
    // printf("here end %d\n",node_size(node_hf));
    free_node(node_hf);
    // //int_rb_tree_free_value(tree_hf);
    int_rb_tree_free(tree_hf);
    int_rb_tree_free(tree_c);
    // printf("hrere free end\n");

    return v2;
}

Mesh* mesh_simplyfied_mesh_boundires(Mesh* m,double threshold,
    int (*vertex_is_delete)(template_v*),
    template_v** vs1,int len)
{
    // printf("begin too simply\n" );
    template_v** vs=(template_v**)malloc(sizeof(template_v*)*len);
    memmove(vs,vs1,sizeof(template_v*)*len);
    // printf("len;%d\n",len);
    for(int i=0;i<len;i++)
    {
        // printf("i :%d ",i);
        if(vs[i]==NULL){continue;}
        // printf("i:%d %d\n",i,vs[i]->id);

        for(int j=0;j<len;j++)
        {
            int k=(j+1+i)%len;
            if(k==i){break;}
            if(vs[k]==NULL)
            {
                continue;
            }
            // printf("k:%d\n",k);

            double dis=distance_of_two_points(vs[i]->point,vs[k]->point,3);
            // printf("dis:%lf k:%d %d\n",dis,k,vs[k]->id);

            if(dis<threshold||
                ( vertex_is_delete!=NULL&&
                    vertex_is_delete(vs[k])==1))
            {
                if(merging_two_vertices(m,vs[k],vs[i])!=NULL)
                {
                // printf("k is null\n");
                    vs[k]=NULL;
                }
            }
            else
            {
                break;
            }


        }
    }
    // printf("end tool simply\n" );

    return m;
}

static double compute_triangle_energy(template_v* v1,
    template_v* v2,
    template_v* v3)
{
    double a= distance_of_two_points(v2->point,v3->point,3),
        b= distance_of_two_points(v1->point,v3->point,3),
        c= distance_of_two_points(v1->point,v2->point,3);

    double cosa = (b*b+c*c-a*a)/(2*b*c);
    double cosb = (a*a+c*c-b*b)/(2*a*c);


    double anglea = acos(cosa) , angleb = acos(cosb);

    return (anglea>=angleb?anglea:angleb);
}

typedef struct Suture_Vertex{
    // 0 是上 ,1是下
    int shang_or_xia;
    int i;
    template_v* v;
}Suture_Vertex;

static inline void suture_vertex_init(Suture_Vertex*sv)
{
    sv->shang_or_xia=0;
    sv->i=0;
    sv->v=NULL;
}

// is_loop： vb1和vb2 是否构成loop
//

static Node* libcell_tools_suture_two_boundarries(
    Suture_Vertex** vb11,int len1,
    Suture_Vertex** vb22,int len2,
    int is_loop)
{
    Node* re=NULL;
    int wei1=0,wei2=0;
    // int len1=0,len2=0;
    if(is_loop)
    {
        double temp_dis=distance_of_two_points(
        vb11[wei1]->v->point,vb22[wei2]->v->point,3);
        for(int i=1;i<len2;i++)
        {
            double temp_dis1=distance_of_two_points(
                vb11[wei1]->v->point,vb22[i]->v->point,3);
            if(temp_dis1<temp_dis)
            {
                temp_dis=temp_dis1;
                wei2=i;
            }
        }
    }

    // printf("len1:%d len2:%d\n",len1,len2);
    int i=0,j=0;
    int thre_len1=len1,thre_len2=len2;
    if(!is_loop)
    {
        thre_len1--;
        thre_len2--;
    }
    while(i<thre_len1||j<thre_len2)
    {
        Suture_Vertex**vs=(Suture_Vertex**)malloc(
            sizeof(Suture_Vertex*)*3);
        // re=node_overlying(re,vs);
        re=node_pushback(re,vs);
        vs[0]=vb11[((wei1+i)%len1)];
        vs[1]=vb22[((wei2+j)%len2)];
        // printf("wei1:%d wei2:%d\n",wei1,wei2);
        if(i<thre_len1)
        {
            if(j<thre_len2)
            {
                double ener1=compute_triangle_energy(
                    vs[0]->v,vs[1]->v,vb11[((wei1+i+1)%len1)]->v);
                double ener2=compute_triangle_energy(
                    vs[0]->v,vs[1]->v,vb22[((wei2+j+1)%len2)]->v);
                if(ener1<ener2)
                {
                    vs[2] = vb11[((wei1+i+1)%len1)];
                    // printf("%d %d vb11:%d\n",vs[0]->id,vs[1]->id,vs[2]->id );
                    i++;
                }
                else
                {
                    vs[2] = vb22[((wei2+j+1)%len2)];
                    // printf("%d %d vb22:%d\n",
                    //     vs[0]->id,vs[1]->id,vs[2]->id );
                    j++;
                }
            }
            else
            {
                vs[2] = vb11[((wei1+i+1)%len1)] ;
                // printf("%d %d vb11:%d\n",vs[0]->id,vs[1]->id,vs[2]->id );
                i++;
            }
        }
        else
        {
            vs[2] = vb22[((wei2+j+1)%len2)];
            // printf("%d %d vb22:%d\n",
            //             vs[0]->id,vs[1]->id,vs[2]->id );
            j++;
        }
        // printf("%d %d %d",)
        // printf("i:%d j:%d\n",i,j);
    }

    return node_reverse(re);
}


static void libcell_fill_vb_index_according_sampling_dis(
    Suture_Vertex**vb,int len,double sampling_dis,
    Suture_Vertex** vb_s,int* lens)
{
    *lens=0;
    if(vb==NULL||len<=0){return;}
    int i=0;
    while(i<len)
    {
        // printf("lens:%d\n",*lens);
        vb_s[*lens]=vb[i];
        *lens=1+*lens;
        i++;
        // printf("1i:%d\n",i);
        for(;i<len;i++)
        {
            double dis=distance_of_two_points(
                vb_s[*lens-1]->v->point,
                vb[i]->v->point,3);
            // printf("dis:%lf\n",dis);
            if(dis >= sampling_dis)
            {
                break;
            }
        }
        // printf("i:%d\n",i);
    }
}


//
//
static Node* libcell_my_get_detail_suture_result(
    Suture_Vertex** vb11,int len1,
    Suture_Vertex** vb22,int len2,
    Node*tris)
{
    Node* re=NULL;
    Suture_Vertex**vs=(Suture_Vertex**)(tris->value);
    int i=vs[0]->i,j=vs[1]->i;
    Node*nit=tris;
    int temp_i=i,temp_j=j;
    Suture_Vertex** vbs1=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len1);
    Suture_Vertex** vbs2=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len2);
    while(nit!=NULL)
    {
        temp_i=i;temp_j=j;
        for(;nit!=NULL;nit=(Node*)(nit->Next))
        {
            // Suture_Vertex**vs=(Suture_Vertex**)(nit->value);
            vs=(Suture_Vertex**)(nit->value);
            if(vs[2]->shang_or_xia==0)
            {
                temp_i=vs[2]->i;
            }
            else
            {
                temp_j=vs[2]->i;
            }
            if(temp_i!=i&&temp_j!=j){break;}
        }
        // printf("i:%d j:%d temp_i:%d temp_j:%d\n",i,j,temp_i,temp_j);

        memset(vbs1,0,sizeof(Suture_Vertex*)*len1);
        memset(vbs2,0,sizeof(Suture_Vertex*)*len2);

        int lens1=0,lens2=0;

        for(int k=0;k<len1;k++)
        {
            int k1=(k+i)%len1;
            vbs1[lens1]=vb11[k1];
            lens1++;

            if(k1==temp_i){break;}
        }

        for(int k=0;k<len2;k++)
        {
            int k1=(k+j)%len2;
            vbs2[lens2]=vb22[k1];
            lens2++;
            if(k1==temp_j){break;}
        }

        Node*noden= libcell_tools_suture_two_boundarries(
            vbs1, lens1,
            vbs2, lens2,
            0);

        // for(Node* nit1=noden;nit1!=NULL;nit1=(Node*)(nit1->Next))
        // {
        //     Suture_Vertex**vs1=(Suture_Vertex**)(nit1->value);
        //     printf("%d %d %d %d %d %d\n",vs1[0]->i,vs1[1]->i,vs1[2]->i,
        //         vs1[0]->v->id,vs1[1]->v->id,vs1[2]->v->id);
        // }


        re=node_plus(re,noden);

        i=temp_i;j=temp_j;
        if(nit==NULL)
        {
            break;
        }
        nit=(Node*)(nit->Next);
    }


    free(vbs1);free(vbs2);
    return re;
}


// 新版本的缝合边界，会对边界进行重采样
// sampling_dis :采样距离
//
Node* libcell_tools_suture_two_boundarriesn(Node* vb1,
    Node* vb2,int is_loop,double sampling_dis)
{
    int len1=0,len2=0;
    len1=node_size(vb1);
    len2=node_size(vb2);

    Suture_Vertex**vb11=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len1);
    Suture_Vertex**vb22=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len2);

    int i=0;
    for(Node* nit=vb1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Suture_Vertex*sv=(Suture_Vertex*)malloc(sizeof(Suture_Vertex));
        suture_vertex_init(sv);
        vb11[i]=sv;
        sv->v=(template_v*)(nit->value);
        sv->i=i;
        i++;
    }

    i=0;
    for(Node* nit=vb2;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Suture_Vertex*sv=(Suture_Vertex*)malloc(sizeof(Suture_Vertex));
        suture_vertex_init(sv);
        vb22[i]=sv;
        sv->v=(template_v*)(nit->value);
        sv->i=i;
        sv->shang_or_xia=1;
        i++;
    }

    // Node* re=NULL;

    Suture_Vertex**vbs1=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len1);
    Suture_Vertex**vbs2=(Suture_Vertex**)malloc(sizeof(Suture_Vertex*)*len2);
    int lens1=0,lens2=0;
    // printf("here1 len1:%d len2:%d\n",len1,len2);
    libcell_fill_vb_index_according_sampling_dis(
        vb11,len1,sampling_dis,
        vbs1,&lens1);
    // printf("here2\n");

    libcell_fill_vb_index_according_sampling_dis(
        vb22,len2,sampling_dis,
        vbs2,&lens2);

    // printf("begin tris\n");
    Node*tris= libcell_tools_suture_two_boundarries(
        vbs1,lens1,
        vbs2,lens2,
        is_loop);


    // printf("end tris\n");
    Node* re=NULL;

    // re=tris;
    re= libcell_my_get_detail_suture_result(
        vb11,len1,
        vb22,len2,
        tris);

    for(Node* nit=re;nit!=NULL;nit=(Node*)(nit->Next))
    {
        Suture_Vertex**vs=(Suture_Vertex**)(nit->value);

        // printf("%d %d %d\n",vs[0]->v->id,vs[1]->v->id,vs[2]->v->id);
        template_v**vs1=(template_v**)malloc(sizeof(template_v*)*3);
        nit->value=vs1;
        vs1[0]=vs[0]->v;vs1[1]=vs[1]->v;vs1[2]=vs[2]->v;

        free(vs);

    }

    free_node_value(tris);
    free_node(tris);

    free(vbs1);free(vbs2);
    for(i=0;i<len1;i++)
    {
        free(vb11[i]);
    }
    for(i=0;i<len2;i++)
    {
        free(vb22[i]);
    }
    free(vb11);free(vb22);
    return re;
}
